wget ftp://ftp.ensembl.org/pub/release-102/gff3/homo_sapiens/Homo_sapiens.GRCh38.102.gff3.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.1.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.2.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.3.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.4.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.5.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.6.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.7.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.8.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.9.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.10.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.11.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.12.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.13.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.14.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.15.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.16.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.17.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.18.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.19.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.20.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.21.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.22.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.X.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.Y.fa.gz
gunzip *gz
cat Homo_sapiens.GRCh38.dna.chromosome.*.fa > hg38.fa
rm Homo_sapiens.GRCh38.dna.chromosome.*.fa
anchorwave gff2seq -r hg38.fa -i Homo_sapiens.GRCh38.102.gff3 -o cds.fa
minimap2 -x splice -t 60 -k 12 -a -p 0.4 -N 20 hg38.fa cds.fa > ref.sam
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/002/915/635/GCA_002915635.2_ASM291563v2/GCA_002915635.2_ASM291563v2_genomic.fna.gz
gunzip GCA_002915635.2_ASM291563v2_genomic.fna.gz
minimap2 -x splice -t 111 -k 12 -a -p 0.4 -N 20 GCA_002915635.2_ASM291563v2_genomic.fna cds.fa > axolotl.sam
perl alignmentToDotplot.pl Homo_sapiens.GRCh38.102.gff3 axolotl.sam > axolotl.tab
anchorwave proali -i Homo_sapiens.GRCh38.102.gff3 -r hg38.fa -a axolotl.sam -as cds.fa -ar ref.sam -s GCA_002915635.2_ASM291563v2_genomic.fna -n axolotl.anchors -R 1 -Q 1 -ns > axolotllog
grep "block begin" axolotl.anchors | wc -l ## 92
changetoM <- function ( position ){
position=position/1000000;
paste(position, "M", sep="")
}
library(ggplot2)
library(compiler)
enableJIT(3)
## [1] 3
library("Cairo")
data =read.table("axolotl.tab")
data=data[which(data$V1 %in% c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10" ,"11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "X", "Y")),]
data=data[which(data$V3 %in% c("CM010927.1", "CM010928.1", "CM010929.1", "CM010930.1", "CM010931.1", "CM010932.1", "CM010933.1", "CM010934.1", "CM010935.1", "CM010936.1", "CM010937.1", "CM010938.1", "CM010939.1", "CM010940.1", "CM010941.1", "CM010942.1", "CM010943.1", "CM010944.1", "CM010945.1", "CM010946.1")),]
data$V1 = factor(data$V1, levels=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10" ,"11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "X", "Y"))
data$V3 = factor(data$V3, levels=c("CM010927.1", "CM010928.1", "CM010929.1", "CM010930.1", "CM010931.1", "CM010932.1", "CM010933.1", "CM010934.1", "CM010935.1", "CM010936.1", "CM010937.1", "CM010938.1", "CM010939.1", "CM010940.1", "CM010941.1", "CM010942.1", "CM010943.1", "CM010944.1", "CM010945.1", "CM010946.1"))
CairoPNG(file="axolotl.png",width = 7200, height = 8000)
ggplot(data=data, aes(x=V4, y=V2))+geom_point(size=4, aes(color=V5))+facet_grid(V1~V3, scales="free", space="free" )+ theme_grey(base_size = 80) +
labs(x="axolotl", y="hg38")+scale_x_continuous(labels=changetoM) + scale_y_continuous(labels=changetoM) +
theme(axis.line = element_blank(),
panel.background = element_blank(),
strip.text.x = element_text(angle = 90),
panel.border = element_rect(fill=NA,color="black", size=0.5, linetype="solid"),
axis.text.y = element_text( colour = "black"),
legend.position='none',
axis.text.x = element_text(angle=300, hjust=0, vjust=1, colour = "black") )
dev.off()
## png
## 2
data =read.table("axolotl.anchors", head=TRUE)
data = data[which(data$refChr %in% c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10" ,"11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "X", "Y")),]
data = data[which(data$queryChr %in% c("CM010927.1", "CM010928.1", "CM010929.1", "CM010930.1", "CM010931.1", "CM010932.1", "CM010933.1", "CM010934.1", "CM010935.1", "CM010936.1", "CM010937.1", "CM010938.1", "CM010939.1", "CM010940.1", "CM010941.1", "CM010942.1", "CM010943.1", "CM010944.1", "CM010945.1", "CM010946.1")),]
data$refChr = factor(data$refChr, levels=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10" ,"11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "X", "Y"))
data$queryChr = factor(data$queryChr, levels=c("CM010927.1", "CM010928.1", "CM010929.1", "CM010930.1", "CM010931.1", "CM010932.1", "CM010933.1", "CM010934.1", "CM010935.1", "CM010936.1", "CM010937.1", "CM010938.1", "CM010939.1", "CM010940.1", "CM010941.1", "CM010942.1", "CM010943.1", "CM010944.1", "CM010945.1", "CM010946.1"))
data$strand = factor(data$strand)
data = data[which(data$gene != "interanchor"),]
p = ggplot(data=data, aes(x=queryStart, y=referenceStart))+geom_point(size=4, aes(color=strand), alpha=0.4)+facet_grid(refChr~queryChr, scales="free", space="free" )+ theme_grey(base_size = 80) +
labs(x="axolotl", y="hg38")+ scale_x_continuous(labels=changetoM) + scale_y_continuous(labels=changetoM) +
theme(axis.line = element_blank(),
panel.background = element_blank(),
strip.text.x = element_text(angle = 90),
panel.border = element_rect(fill=NA,color="black", size=0.5, linetype="solid"),
axis.text.y = element_text( colour = "black"),
legend.position='none',
axis.text.x = element_text(angle=300, hjust=0, vjust=1, colour = "black") )
CairoPDF(file="axolotl.anchors.pdf",width = 90, height = 100)
p
dev.off()
## png
## 2
CairoPNG(file="axolotl.anchors.png",width = 7200, height = 8000)
p
